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Abstract 

In this paper we find an analytical solution of the equilibrium ion distribution 
for a toroidal model of a ionic channel, using the Perfect Screening Theorem 
(PST)[[]. The ions are charged hard spheres, and are treated using a scaling 
Mean Spherical Approximation (SMSA) 0. 

Understanding ion channels is still a very open problem, because of the 
many exquisite tuning details of real life channels. It is clear that the electric 
field plays a major role in the channel behavior, and for that reason there has 
been a lot of work on simple models that are able to provide workable theories. 
Recently a number of interesting papers 0] |SJ El HI EH EH have appeared 
that discuss models in which the effect of the geometry, excluded volume and 
non-linear behaviour is considered. 

We present here a 3D model of ionic channels which consists of a charged, 
dcformablc torus with a circular or elliptical cross section, which can be flat or 
vertical (close to a cylinder). Extensive comparisons to MC simulations were 
performed. 

The new solution, which is simple and accurate, opens new possibilities, such 
as studying flexible pores |12| and water turning phase transformations inside 
the pores, using an approach similar to that used on flat crystal surfaces [T^1H*4] 

1 Introduction 

We dedicate this contribution to Prof. Ben Widom, one of the leading figures 
in Statistical Mechanics. 

The study of the transport process in membrane ion channels is complicated 
by the presence of the protein walls, the interaction with ions, water molecules 
and the electric field profile, which determine many of the salient properties of 
ion channels [S]- Computing the electric potential profile everywhere in a real 
channel is difficult, if not impossible, because of the complexity of the system, 
and for that reason simplified models have been used: Kuyucak et al 0] 03 , have 
studied circular toroidal channels using the linear Poisson-Boltzmann equation, 
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in which the ions are treated as point charges. Excluded volume effects have 
been included [5] to explain ion selectivity. Furthermore nonlinear effects, which 
are important, are included in the ID, non-linear Poisson Boltzmann PNP mod- 
els of Eisenbcrg ct al.|SlE3El- Excluded volume effects come into play when 
molecular solvents are used [51 [7|. Recently the effects on porin arising from the 
rotation of water molecules were discussed [T^]. 

We propose here a SMSA solution |5j of a ion channel model which consists 
of a toroidal ring with either circular or elliptical (prolate or oblate) cross sec- 
tion. The solution is given in terms of a few MSA scaling parameters provide an 
accurate description of the charge disribution and can be obtained from a vari- 
ational theory. The major feature of our solution is that it satisfies the Perfect 
Screening Theorems (PST)pQ. This is not only a physical requirement, but also 
a technical advantage, because it can used to include discrete molecular solvents 
such as hard dipoles ^H], an d water ^l- This has been used in the theory of a 
phase transition that involves the turning of water by an electric field |14j. 

The non- linear Poisson-Boltzmann case has been discussed elsewhere |17II18| 
and as a matter of fact, is implicit in our present work. 

To study the dynamics of ions in a channel, one needs to compute the forces 
acting on each of the ions , including mobile, induced and fixed charges and 
the applied electrical field. This could be coupled with Brownian dynamics 
simulations or other coarse grained simulations. Because this computation has 
to be repeated at every step, the existence of analytical solutions in a relevant 
geometry is imperative for simulations at realistic time scales. 

1.1 The perfect screening sum rule (PST) 

One remarkable property of mixtures of classical charged particles is that be- 
cause of the very long range of the electrostatic forces, they must create a 
neutralizing atmosphere of countcrions, which shields perfectly any charge or 
fixed charge distribution. Otherwise the partition function, and therefore all 
the thermodynamic functions, will be divergent £Q. This sum rule is intuitive 
and widely accepted for spherical systems. But it is also true for non-spherical 
systems. It has been explicitly verified for non spherical systems in simulations 
and also using exact results for the Jancovici model |19l I2(J| . 

For spherical ions this means that the internal energy E of the ions is always 
the sum of the energies of spherical capacitors. For any approximation the exact 
form of the energy is 



j3 = 1/kT is the usual Boltzmann thermal factor, £ is the dielectric constant, e 
is the elementary charge, and ions i have charge, diameter and density Zi, a%, 
pi, respectively. Ti is the shielding length for ion i. 
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An equivalent electrostatic model of eq.© is the ion immersed in a con- 
ducting media, which means that the energy will depend only on the charge 
of the ion and not on its internal distribution: The energy will be the same if 
the charge is concentrated at the center of the sphere or uniformly distributed 
over the surface of the ion [^[221 [23]: Clearly the potential at the surface has 
to be constant and the proper boundary conditions for the electrostatics are 
the Dirichlct boundary conditions. The same considerations apply for an arbi- 
trary fixed charge distribution. In the case of a solid toroid of circular section 
immersed in a conducting media the potential outside of the torus will be the 
same if we took a circular wire at the center of the torus or a charged ring with 
a constant surface potential. We remark that this is true strictly for Dirichlet 
boundary conditions. The more general case is more complicated. 

For non spherical systems the perfect screening theorem requires that all 
multipoles of the fixed charge distributions be compensated by the mobile charge 
distribution. As we will sec below this implies a very substantial simplification 
of the solution of linear PB equation since every multipole of the countercharges 
distribution cancels the fixed charges multipoles and all cross terms are zero. 



2 The charged torus 

Analytical solutions of closures of the Orstein-Zernike (OZ) or Wertheim-Ornstein- 
Zernike (WOZ) equations are only possible in odd parity spaces. The torus is an 
even parity, 2-dimensional object, and for that reason there is no direct analyti- 
cal solution possible of the MSA or LPBE. The expansion in spherical harmonics 
on the other hand is always possible and for our model is convergent within a 
reasonable (even small!) number of spherical harmonics j^U- The beauty of the 
PST is that it de-couples all the multipole terms, and for that reason we are 
able to solve the LPBE to all orders in closed form. 

Consider first a circular section b = b rea i + a/2 torus of radius a. The 
diameter of the ions is a. The electrostatic equivalent system is a ring wire of 
radius d = a and the torus immersed in a conducting media. Poisson's equation 
for the potential in a charged system is |2U 

V 2 ^(r) = %(r), (2) 

Here the charge density q(r) at r is the sum of the fixed ring and the mobile ion 
charges 

q(r) = qrvng 0) + ^ & ( r ) ( 3 ) 

i 

where </>(r) is the potential at r = R, z . 
The formal solution of this equation is 
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The potential generated by a charge Q on the ring of radius d is given by 

2Q , . ARd 

^ rmg (R,z)^ e[{R _ d)2 + z2] K(- m ); rn= {R d)2 + z2 (5) 

where K(m) is the elliptic function 

K{m) = [* d \ , (6) 

Jo v 1 — m z sin 

which satisfies the homogeneous Poisson equation 

V 2 = 0. (7) 

However, the inhomogeneous Poisson equation @ has no closed form analytical 
solution. We need to expand our problem in a suitable basis. We use spherical 
harmonics expansion because its good analytical behavior (it has been exten- 
sively used in astrophysics but more important, because of the PST P 
the terms are decoupled to each order in the expansion. For the ring source po- 
tential there are three regions which correspond to outside and inside a sphere 
of radius a of the ring. 

Mr) = <P ext e Heav {r-a-b) + 4> rin ^ Heav {a + b-r)e Heav {r-a + b) 
+ 4> mt 6 Heav (a-b-r) 

(8) 

The parameter b is the effective diameter of the circular torus. We assume that 
the ions are of diameter a and therefore the real radius of the toroid is 

breal — b — — 

When the width of the ring is b = we get 

oo 

(j) ext {r) =^2P e (cos9)r~ {e+1) M^ t (9) 



with 



and 



with 



M ext = !±p t (o)a*. (io) 

OO 

nt (r) = r e Pi(cos6)Mi nt (11) 

M int = Qp e ( ) a -(z+i). (12) 
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The fixed multipole moments Me are the same for a ring of charge Q or a solid 
toroid with uniform potential on its surface ( Dirichlct boundary conditions). 
This will be true for a pore immersed in a conducting electrolyte. We will take 
advantage of this fact computing the moments for the ring of zero with, and 
then using them in the calculation of the potential for a toroid of finite width. 

The calculation of the multipole moments Me for the general case of ellipti- 
cal toroids is left for a future publication. 





Fig. la. Charged rings with elliptical cross sections. 
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Fig.lb. 



3 The Solution of the Linear Poisson Boltzmann 
Equation (LPBE) and the Scaling MSA 

But the LPBE is accurate only for very dilute solutions because the ions are 
point charges. The MSA is the LPBE, but with the mathematically correct 
treatment of the excluded volume effects. The MSA which is derived from the 
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Mean Spherical Model ^ provides a coherent and simple description 

of the properties of a large class of systems in terms of a very small set of 
scaling parameters T a This includes ionic solutions, water and polyelectrolyte 
|2IlE21ESIElESlEniEZlEHllSni- The functional form of the thermodynamic 
parameters in the different modified MSA theories is the same as in the LPBE, 
however the excluded volume is treated exactly in the MSA which satisfy the 
Onsager high density bounds 00] 123 113 i and is asymptotically exact at large 
concentrations. One can show that simple transformations lead to the proper 
high density behavior. The Debye screening length k in the DH theory becomes 
the MSA screening length T 



i 



^£"»=j =>■ r4*«»)-i (13) 



It can be shown that the proper high density behavior in the MSA stems 
from the fact that the entropy is of the form 

AS (MSA) = _ ks ^_ (M) 
37T 

where ks is the Boltzmann constant. For nonsphcrical systems this generalizes 

tong 

AS (MSA) = _ kB J2^C (15) 

x=-t 7F 

where x is the index of the irreducible representation |41j and £ is the order of 
the spherical harmonic in eq. © . This immediately suggests 03 13 that T x can 
be determined by the variational expression 

d[(3AA x (T x )} _ d[(3AE x (T x ) + T 3 x /(37T)} 

dr x dr x { ] 

For the simple restricted case of an equal size ionic mixture we get equation 
qi3|). For more complex systems, like the general polyelectrolyte this equation 
is a new relation to be found. For flexible polyelectrolytes it has been derived 
from the binding MSA (BIMSA) [321131150] . 

In this work we use the LPBE in conjunction with the perfect screening 
theorem (PST)PP, to derive the functional form of the solution to our problem. 
This functional form provides an astonishingly simple and good representation 
of our simulation data. From cq. (|16f) we get a very good first approximation to 
the charge distribution obtained from our extensive simulations. 
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3.1 Solution of the LPBE 

Consider eq. Q : In the linear Poisson Boltzmann approximation we write as a 
convolution: 

k 2 (Y 



W=0o-^ ( -) *0(r) (17) 



Taking the Fourier transform of both sides 



K 2 (\ 



where 



Then. 



where 



^( fc )-^-^(-)0(fc) (18) 



A- 



2 



m = pt^ 00 (20) 



00 = <Po + +00 



7^0 

Using Raylcigh's formula 

e lkrcose = J2(2£+l)i e P e (cos6)j e (k 



we get 



rin t 



and 

J+2 



tint 
b 



1=0 

Lext 



(21) 



*«-FT^* + PT^* , + FT^ (22) 



= 4Tr^ Mi nt P e {cos9)i e j e (kr)r e+2 dr, (23) 



4^5]^ Pt(cos6) i l ^- i/+i(fto). (24) 
We find 0q x4 from equation (13): 

<^(fc) = 4tt^ Mr* P/(co*0) l f ^r- j/+i(fco). (25) 

£=0 

and 0q 1 " 9 is a linear combination of <pQ Xt and 0™' with coefficients to be deter- 
mined by the boundary conditions. 

The inverse Fourier transform of 4> ext (r) is 

4> ext {r) = 4^P,(c OS 0)i?f '(r), (26) 
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where, after some calculations [52] 



9 f°° k 3 

Mr~a t+2 / dk-f-^ U {kr) H+1 {ka) 



= Ml nt n 2 a e+2 i t {nr) k l+1 {na) 

(27) 

Similarly 

= E^v-v^^ d 3 fc __ J£ _ i(fca) ^ 



ik.r 



(28) 



Finally we find the LPBE potential 

<t> ex \r) = E^( cos0 )[^ k M\ 

i 

<f> int (r) = J2 P dcos9)[Aeh(Kr)} 

1 

f iB9 (r) = ^Pi(co89){Cti t (Kr)+Dtkt(Kr)} 

£ 

where 



(29) 



At = Mj nt n 2 a e+2 k i+1 ( K a), (30) 

B e = M e t xt k 2 a 1 ' 1 i e+1 (na) (31) 

and Ct and Dg are found from the boundary conditions on the spheres of radius 
r = a±b. The potential at the interface is constant, so that in the simplest case 
they are obtained by matching the normal derivatives at the toroid's surface. 
This means that the matching at the toroid's surface must include a geometrical 
factor 



b(z) = \Jb 2 - z 2 ; b> z 



z 



cos 9 = (32) 

yja 2 + b 2 ± 2ab{z) 

This factor only matters when z m b. 

From the simulations we know these coefficients depend also on nonlinear 
effects. But what is clear is that they provide a good representation of the re- 
sults of our current simulations. 



The spherical bessel functions ii(nr) and kg(nr) satisfy the proper 
boundary conditions : 

<hnt(r, 6) = A l Pi{cos9) i<(«r), (0 < r < a - b) (33) 
i 

The exterior potential is 

<l>ext (r, 9)=Y,B e P t {cos6) ki{nr), (r>a + b) (34) 
l 

As was discussed elsewhere El 122 it has been shown that in the many cases 
where an analytical solution of the MSA, for complex systems is available, the 
solution can be obtained from a simple variational principle. The actual solution 
of the MSA is more complex for higher values of I. But in our simulations the 
largest term is the one with I = so that ea. (|13J) is a quite decent approximation. 
Then our expression is 

MR, z) = <t> ext Hea v(r -a- b(z)) + <f> int 6 Heav (a - b(z) - r) 
+<f> rin9 e Heav {a + b{z) - r)6 Heav {r -a + b(zj) 

(35) 

with 

<P ext (r) = ^2P e (cos0)[B t kt(Tr)] 

i 

r'(r) = ^Piicose^Ae u(Tr)] 

t 

0™9(r) = Y j P l {cos6){C l u{Yr)+D l k t {Yr)} 

i 

where 

At = AirMi nt r 2 a t+2 k e+1 (Ta), (36) 
B t = AttM^T 2 a 1 " 1 i t+1 {Ta) (37) 

4 Computer Simulations 

The geometry of the system is displayed in Figure 1. We use a cylindrical ring 
of radius 6=1/2 The size of the ions is also a — 1. We have chosen this model 
because it is closest to actual channel geometry [21 0] . As is shown in figure 1 
we can stretch or flatten the torus so that it can be a narrow channel or just 
a pore in a flat membrane. In our solution this means that we merely have to 
increase the number of multipoles in the expansion eq.©. As was shown by 
Tohlinc et al. the expansion of toroidal systems using spherical harmonics can 
be carried out up to quite large values of i ~ 1000 |2SI - Simulations discussed 
in this paper are only for cylindrical toroids. We simulate a cylindrical section 
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torus for which only a few terms are needed (typically 3 or 4). We use standard 
Monte Carlo simulation techniques [31 023 El . Our entire system is confined to 
a large sphere, and our charged ring is lying flat in the x-y plane, so that z = is 
the plane of the ring. Typically the radius of the sphere is 20-26 times the ionic 
diameter, although runs have been performed with larger spheres to assure that 
the influence of the walls is small. The number of ions varied between 200 and 
512. The sampling bins were cylindrical sections, which makes the data near 
the axis of the torus particularly noisy. For that reason we have used extremely 
long runs ( 10 8 steps in some cases) and even so we had to discard the noisy 
central region. Independence of the size of the spherical container was verified. 

In the simulations the ring had a radius of 3, and since the ions had a radius 
of 1, the effective exclusion region goes from 2 to 4 units. The results of figures 
2 and 3 for the charge distributions in the x-y plane are well represented by 
eq.JSHJ), using the first two terms in the expansion. Figures 4 and 5 represent 
the same charge distributions for different values of z. The curve for z=l marks 
the transition region where the ring ends. The curve is a perfect parabola given 
by eq. (32) and spans the region a ± b(z). For z > b = 1 the profile is simply 
a superposition of the interior and exterior regions. The overall agreement of 
these figures with the theory is comparable to those of figures 2 and 3. 

Figures 6 and 7, which correspond to z = 0, and charges Q = 4, 12, 20, are 
also in excellent agreement with the theory of eq. (34), in the sense that the 
functional form is always the same. However the parameters Ai, Bi,Ci, Dg , 
have a nonlinear dependence on the charge Q. We hope to discuss this point in 
the near future. 

Our simulations are designed to satisfy the perfect screening theorem, £Q, 
and therefore we will not have large multipoles (I > 2) induced by the periodic 
boundary conditions. 
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Fig. 2. Comparison of theory and computer simulation of the anions density 
for a ring of internal diameter 3, charge 20. The dotted line is the theory 
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Fig. 3. Comparison of theory and computer simulation of the cations density 
for a ring of internal diameter 3, charge 20. The dotted line is the theory. 
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Fig. 5. g a nion(R, z) for the same parameters as in figure 4. The graph also 
shows three regions of ea. (|21|l when z < b and the two regions for z > b 
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Fig. 6. Charge dependence of g a nion(R, z) for z=0. 
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Fig. 7. Charge dependence of g ca tion(R, z) for z=0. 
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5 Discussion of Results 



This research is part of an ongoing effort to describe complex systems by means 
of a small set of scaling parameters. The basic idea is to use a singlet density de- 
scription combined with exact relations, such as the Onsager asymptotic limits, 
or the Wertheim association limits or, as in the present case, the perfect screen- 
ing theorem. For example we have developed a theory of flexible polyelectrolyte 
of arbitrary length that satisfies an exact relation for infinite length [T2l 1 18| . and 
for that reason agrees extremely well with computer simulations. The same 
theory has been used to describe a completely flexible ring (Bernard and Blum, 
unpublished): In this case the 'channel' is a chain of spherical beads of arbitrary 
charge and diameter. The toroidal channel is a model that in a certain sense 
represents the other extreme situation since it is rigid and structure-less. The 
'true' model lies somewhere in between. 

Another situation where our solution could be used is to describe the case 
when the water molecule turns, as was discussed by Tajkhorshid et al. |13| . In 
fact we have used our methods to describe a phase transition which is actually 
due to the turning of water in an external electric field |14| . 

Also dielectric effects can be included, adding another scaling parameter that 
has the meaning of a polarizability. This was recently published ^5] . Actually 
equation (16) is derived in that paper. And finally, external fields are easy to 
include and the effect of membranes can be mimicked by a suitable deformation 
of the torus shown in figure 1. 

We have presented here a theory of the equilibrium distribution of ions in 
a charged channel that is surprisingly simple. The fact that so few parameters 
can describe the distribution of charges inside and outside the channel is very 
appealing. A feature of our approach is that it can be easily extended to include 
polarization effects, solvent and discrete structures. Because of the availability 
of an analytical solution of the octupolar model of water ^1] it is another defi- 
nite possibility. 
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